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ABSTRACT 

An extremely fast time-harmonic finite element solver developed for the transmission analysis of photonic crystals 
was applied to mask simulation problems. The applicability was proven by examining a set of typical problems 
and by a benchmarking against two established methods (FDTD and a differential method) and an analytical 
example. The new finite element approach was up to 100 x faster than the competing approaches for moderate 
target accuracies, and it was the only method which allowed to reach high target accuracies. 
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1. INTRODUCTION 

The complexity of modern photolithography makes extensive simulations indispensable.^ Modern lithography 
simulators include modules describing illumination, transfer of the optical field through the mask and aberrating 
•i^ ' optical system of the lithographic equipment, the propagation inside the photoresist, the processes leading to 
the resist image and - in advanced systems - the etching processes leading to the etched image. After nearly 
^ two decades of lithography simulation, most of the modules along the simulation chain have attained a high 
■ " " ' degree of maturity. However, the simulation of light propagation through phase shift masks, also applied for the 
stand-alone analysis of masks and mask tolerances, is still challenging in terms of both computational time and 
accuracy of the results. 

The computation of the print image of a whole chip remains extremely demanding although approximations, 
multi-threading and even hardware accelerators are applied to reduce the runtime of simulations. Rigorous 
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simulations arc restricted today to small areas and even those simulations suffer from the high computational 
effort. At the same time, the progress on the semiconductor roadmap forces the need of rigorous 3D simulations, in 
particular also for alternating and attenuated phase masks. Experimental investigations of the polarization effects 
in Hyper NA immersion lithography^ support this assertion. Further, the demand to assess the process stability 
by exploring several dose/defocus conditions in the process window sharpens the shortage of computational 
resources. 

Keeping this background in mind, we evaluated a frequency- domain finite-element method (FEM) solver for 

Maxwell's eqiiations which has been successfully applied to a wide range of electromagnetic field computations 
including optical and microwave waveguide structures, surface plasmons, and nano-structured materials. In 
addition, the activity was motivated by a preceding, successful benchmarking^ against MPB, a widely used and 
highly sophisticated plane-wave solver used for the Bloch mode analysis of photonic crystals. 

In this contribution, the new FEM solver is benchmarked against an analytical result which is fairly realistic 
for lithography, and against two competing algorithms commonly applied for the simulation of "thick" phase 
masks. For the simulation of periodic mask patterns in lithography the most prominent rigorous simulation 
methods include the finite-difference time domain algorithm (FDTD)^'^ and the modal methods such as the 
differential method^' ^ or the closely related rigorous coupled wave analysis (RCWA).^ The methods differ in 
the way Maxwell's equations are numerically solved and how the boundary conditions of the interfaces to the 
unbound regions above and below the mask are established. We will give a brief description of the algorithms to 
give a rough idea about the relevant parameters influencing simulation speed and accuracy in the next paragraph. 
A brief description of FEM will be given in Section 2. 

The FDTD approach discretizes Maxwell's equations in both time and space and solves the scattering prob- 
lem by simulating the field evolution through time until the time-harmonic steady-state solution is reached. The 
interfaces to the unbound regions are formed by perfectly matched layers (PML). Space and time discretization 
are interdependent ("magic steps"). Speed and accuracy of the simulation depend on the space and time dis- 
cretization, the total time period and on the PML parameters. The differential method describes the propagating 
fields inside the mask by a plane wave expansion. Maxwell's equations are thus converted into a system of linear 
ordinary differential equations (ODEs) relating the scattered waves at the upper and lower mask interfaces. 
Speed and accuracy depend on the number of plane waves and on the resolution used for ODE integration. 

This paper is organized as follows: Section 2 introduces the concept of the FEM solver. Section 3 presents 
several proofs of applicability including 3D simulations, 2D simulations for light scattering off line masks under 
conical incidence, and adaptive refinement of FEM meshes. Section 4 shows benchmark results of the three 
different solvers for a mask with dense lines and spaces illuminated at normal incidence. Section 5 verifies the 
accuracy of the FEM solver by examining a closely related toy example offering an anlytical solution. 



2. FREQUENCY-DOMAIN FEM ANALYSIS FOR PHOTOMASK SIMULATIONS 

This paper considers light scattering off a system which is periodic in the x— and y— directions and is enclosed 
by homogeneous substrate (at Zsub) and superstrate (at Zsup) which are infinite in the — , resp. -|- 2;— direction. 
Light propagation in the investigated system is governed by Maxwell's equations where vanishing densities of free 
charges and currents are assumed. The dielectric coefficient e{x) and the permeability iJ,{x) of the considered 
photomasks are periodic and complex, e{x) = £{x + a), n{x) = ix{x + a). Here a is any elementary vector 
of the periodic lattice. "'^'^ For given primitive lattice vectors ai and 02 an elementary cell C M'^ is defined 
as r2 = € I a: = aiai -|- 0202; < ai, 02 < l} x [Zgub, Zsup\- A time-harmonic ansatz with frequency to and 
magnetic field H{x,t) = e~"^*iJ(x) leads to the following equations for H{x): 

• The wave equation and the divergence condition for the magnetic field: 

V X V X H{x) - J^ii{x)H{x) = 0, f e O, (1) 
e{x) 

W-fj.{x)H{x) = 0, x&n, (2) 



• Transparent boundary conditions at the boundaries to the substrate (at Zsub) and supcrstratc (at Zgup), 
dQ, where -ff" is the incident magnetic field (plane wave in this case), and n is the normal vector on dfl: 

^ V X (i? - i?'") ) X n = DtN{H - i?*"), x G 90. (3) 



£(f) 

The DtN operator (Dirichlet-to-Neumann) is realized with the PML method. This is a generalized 
formulation of Sommerfeld's radiation condition; it can be realized alternatively by the Pole condition 
method. 

• Periodic boundary conditions for the transverse boundaries, dQ, governed by Bloch's theorem^°: 

H{x) = e*^'^w(f), u{x) = u{x + a), (4) 
where the Bloch wavevector G is defined by the incoming plane wave Jf*". 

Similar equations are found for the electric field E{x,t) = e~^'^*E{x); these are treated accordingly. 

The finite-element method solves Eqs. (1) - (4) in their weak form, i.e., in an integral representation. The 
computational domain is discretized with triangular (2D) or tetrahedral (3D) patches. The functional spaces 
are discretized using Nedelec's edge elements, which are vectorial functions of polynomial order (typically second 
order) defined on the triangular or tetrahedral patches. ^'^ In a nutshell, FEM can be explained as expand- 
ing the field corresponding to the exact solution of Equation (1) in the basis given by these elements. This 
leads to a large sparse matrix equation (algebraic problem). For details on the weak formulation, the choice 
of Bloch-periodic functional spaces, the FEM discretization, and our implementation of the PML method we 
refer to previous works. To solve the algebraic problem on a personal computer either standard linear 
algebra decomposition techniques (LU-factorization, e.g., package PARDISO^^) or iterative methods^'' are used, 
depending on problem size. Due; to the use of multi-grid algorithms, the computational time and the memory 
requirements grow linearly with the number of unknowns. 

From the users's point of view, the FEM approach presented here offers the following advantages: 

• The expansion into localized functions (shared with any FEM and FD approach) is adequate for step index 
profiles occuring in masks. 

• The flexibiltiy of triangulations (shared with any FEM approach) allows for the simulation of mask imperfec- 
tions such as sloped etch profiles and for adaptive mesh-refinement strategies leading to faster convergence. 

• The frequency domain approach (shared with any PW method) is adequate for monochromatic or nearly 
monochromatic illumination. 

• Edge elements provide "built-in" dielectric boundary conditions crucial for a high precision simulation of 

step index profiles. 

• The mathematical structure of the algebraic problem allows for the use of very efficient numerical solvers, 
i.e., numerical methods where the computational effort grows linearly with the number of unknowns only. 

• The FEM discretization is characterized by two parameters, the mesh width h and the thickness of the 

PML layer p. It is mathematically proven that the FEM approach converges with a fixed convergence rate 
towards the exact solution of Maxwell-type problems for mesh width /i ^ 0, and p oo}^' This allows 
to easily check whether the attained results can be trusted. 

These advantages result not only in an increased attainable accuracy, but - via the reduced number of unknowns 
- also in a significantly reduced computational effort at moderate target accuracies required for lithography 
simulation. The investigated FEM solver JCMharmony includes adaptive grid refinement, higher order, 3D 
Nedelec elements, advanced periodic and transparent boundary conditions and flexible interfaces to the drivers 
and for postprocessing. Typical computation times for 3D problems are 30 seconds for problems with A'' f« 30 000 
unknowns and 5 minutes for problems with N « 150 000 unknowns solved on an actual standard 64 bit personal 
computer {AMD Opteron). Typical computation times for 2D problems are given in Table 2. 



3. FEATURES OF THE FEM SOLVER 



The range of applications of the FEM approach was examined by means of several characteric tasks in mask 
simulation. 

3.1. Conical Incidence 

Here we investigate line masks illuminated under conical indicence (i.e., oblique incidence with respect to both 
mask plane and grating lines). This is crucial for the accurate simulation of off-axis source points for dipole, 
quadrupole or annular illumination. Figure 1 shows the schematics the geometry of the problem. The geometrical, 
material and source parameters are given in Table 1 (data set 1). The geometry does not depend on the y- 
component, therefore Eqn. (1) reduces to a simpler equation where 2D differential operators act on the 3D 
electric, resp. magnetic fields. Nevertheless, the problem is simulated without any approximations. 
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Figure 1. Schematics of the geometry of a periodic linemask: The computational domain consists of a line of width w 
(at center of the line), height h and sidewall angle /3, on a substrate material Si02, surrounded by air. The geometry is 
periodic in a;-direction with a pitch of px and it is independent on the j/-coordinate. The refractive indices of the different 
present materials are denoted by ni (line), n2 (substrate) and ns (air), ns = 1.0. 
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Figure 2. Diffraction of a plane wave under conical incidence. Intensities of the transmission to the central diffraction 
orders (T = |A(fci)p,i = — 1,0, +1) in dependence on the inclination angle a, where A(k^ is defined in Eqn. (5). The 

angle of rotation is constantly B — 20°. 

The performance of the FEM solver is demonstrated by a parameter scan for varied inclination angle a of 
the source (S'-polarization). The wavevector k of the incident plane wave is attained by a rotation of the vector 
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Table 1. Parameter settings for the simulations in Sections 3.1 (data set 1), 3.2 (set 2), 3.3 (set 3) and 4 (set 4). 



(0,0,27r/A) (where A = Xo/n2, vacuum wavelength Aq, refractive index 712) around the y-axis by the inclination 
angle a and a subsequent rotation around the ^;-axis by the rotation angle 6. In this scan we fix the rotation, 
G = 20° and vary the inclination, a = —25 . . . 25, further parameters are given in Tabic 1. This yields an incident 
wave vector which is scanned from k « (-2.021, -0.736, 4.612)107m to fc « (2.021, 0.736, 4.612)107m. Figure 2 
shows the normalized magnitude of the Fourier coefficients corresponding to the zero and first diffraction orders 
of the scattered light field in dependence on the angle of incidence. A typical computation time for a single 
data point in this scan was 15 sec {N 3 x 10^ unknowns, adaptive grid refinement, computation on a personal 
computer /Intel Pentium IV, 2.5 GHz), resulting in a total time of roughly Ih for the scan with 200 data points. 

3.2. Degree of Polarization of Light Transmitted through a Line Mask 

We have performed a scan over different geometrical parameters by varying the pitch and the linewidth of a 
line- mask. Geometrical, material and source parameters are again given in Fig. 1 and in Table 1 (data set 2). 
Since the plane of incidence is normal to the grating lines {ky = 0), TE- and TM-polarization is supported, i.e., 
the problem becomes scalar. 
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Figure 3. (a) Degree of polarization in dependence on mask pitch {px) for light diffracted to the zero (i = 0) and first 
(i = 1) diffraction order. Several strong Wood's anomalies are indicated at Px = NX. (b) Enlargement of a detail. 



Figure 3 shows the degree of polarization of the zero and first transmitted diffraction orders, defined as 
DOP = {Ite — Itm)/{Ite + Itm) in dependence on the pitch px- The strong Wood anomalies,'' some of 
which are also verified by Tcniber ct al.,^ are caused by the excitation of waves traveling along the mask surface. 
We have constructed transparent boundary conditions for the whole range of investigated pitches, including 
'regular' regions where the transmitted diffraction orders correspond to plane waves with a nonzero ^;-component 



of the wavcvcctor and regions close to A^^ood's anomalies where eertain diffraction orders cannot propagate as 
plane waves anymore. We are currently implementing an adaptive strategy for the PML implementation of the 
transparent boundary conditions in order to automatically account for such effects. 

The average computation time for a single data point in this scan was about 3.5 sec {N ^ 1 — 4 x 10* unknowns, 
depending on geometry size, yielding a relative error of the diffraction intensities of less than 1%; computation 
on a personal computer /Intel Pentium IV, 2.5 GHz), resulting in a total time of roughly 75 min for the scan with 
1200 data points. Similar results are obtained with the solvers SOLID E and Delight (see Section 4), however, 
the Wood's anomalies are not or less accurately resolved. 

3.3. Adaptive Grid Refinement 

By refining the resolution of the geometry-triangulation the accuracy of the solution is increased. FEM-solvers use 
as a standard a regular grid refinement, i.e., in 3D each tetrahedron of the discretization is refined to eight smaller 
tetrahedra, in 2D each triangle of the discretization is refined to four smaller triangles. However, FEM meshes 
also allow for adaptive strategies where only certain elements of the triangulation are refined. The investigated 
solver JCMharmony uses a residuum-based error-estimator^^ for adaptive grid refinement. Obviously, adaptive 
grid refinement is especially useful when the sought solutions are geometrically localized, or when the geometry 
exhibits sharp features, like discontinuities in the refractive index distribution. 
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Figure 4. Comparison adaptive vs. regular FEM grid refinement. Relative error of the transmission intensity in the zero 
(circles) and first (squares) order of TM-polarized light incident onto a periodic line mask, calculated with the FEM solver 
JCMharmony. Solid lines correspond to automatic adaptive refinement, dotted lines correspond to regular refinement. 

We compare the different refinement strategies by observing the convergence of the solutions obtained for 
increased grid resolution (i.e., increased number of imknowns). The setting is similar to the one in Section 3.2. the 
parameters are given in Table 1 (data set 3), and the incident light is TM-polarized. Fig. 4 shows the convergence 
of the relative error of the light intensity, A7 = |/jv,i — Iin{,i\/Iin(,i, in two different diffraction orders for adaptive 
and for regular grid refinement. Here, Ina denotes the light intensity in the diffraction order calculated from 
a solution with N unknowns, /inf,i denotes the intensity calculated on a very fine grid {quasi-exact solution). 
In this example, the use of the error estimator and adaptive refinement yields two orders of magnitude in the 
accuracy of the error for a number of unknowns oi N ^ 10^. 

3.4. Fully 3D Simulations 

As a true 3D example, we have examined the transmission through a mask with a periodic 3D pattern as shown 
in Figure 6a ( "chequerboard pattern"). It is illuminated by a plane wave incident from top. The parameters 
h, ni, 712, ri3, and Aq are the same as in the previous examples. The mask is discretized by a tetrahedral 
mesh supporting second order Nedelec elements. The inital triangulation is shown in Figure 6b. After two 





Figure 5. Unit cell (920 nm x 800 nm x 85.4 nm) of a periodic mask pattern (a) built up from absorbing material (dark 
gray, ni) and air (light gray) and its initial tctraliedral triangulation (b). The mask is illuminated at normal incidence. 
The graphs (c) and (d) show the (projected) vectorial solution and the corresponding intensity gray scale map on a cross 
section well below the output side of the mask, (e) and (f) show intensity maps for cross sections in the middle and at 
the input side of the mask, respectively (white: low intensity, black: high intensity). (See original publication for images 
with higher resolution.) 

refinement steps, the discretization led to a linear system with 2.7 • 10^ unknows. A cross section of tlie 3D 
vectorial solution (c) and cross sections through the 3D intensity distribution (d-f) are also shown in Fig. 5. One 
cleary observes the expected discontinuous behavior of the electric field at material interfaces. For more details 
on 3D FEM computations we refer to a previous work.^° 

4. BENCHMARK OF DIFFERENT RIGOROUS METHODS 

We have performed a benchmark of the previously described FEM solver JCMharmony and two other advanced 
methods, which are also commercially available: The Finite-Difference Time-Domain solver SOLID E,^^ and the 
solver Delight,* which relics on a waveguide-method. In this benchmark we have investigated light propagation 
through a phase mask. The geometry of the problem is outlined schematically in Figure 1. A plane electromag- 
netic wave is incident onto the computational domain with a wavelength of A, a wavevector of = (0, 0, -|-27r/A) 
and a polarization of ^ = {0,Hy,0), with Hy = rul'^ ■ 1 A/m (TM), resp. E = {0,Ey,0), Ey = n^^''^ ■ IV/m 
(TE). The geometrical and material parameters are denoted in Table 1 (data set 4). 

For a quantitative assessment of the different simulation methods we investigate the internal convergence 
behavior of the simulation results, i.e., the deviation of the simulation results from the highest effort result 
obtained with the same method. It turncs out that in all cases the accuracy of the field representation inside the 
computational domain dominates the convergence speed. For FDTD both space and time resolution are affected 
since the stability of the method requires aligning time resolution with spatial resolution. For the differential 
method, the spatial resolution is determined by the number of Fourier coefficients - since the grating considered 
here exhibits discontinuities in the material distribution, evanescent waves play an important role for correctly 
approximating the fields in the mask. For the FEM solver, the accuracy of the solution depends on the number 
of unknowns which is given by the number of geometrical patches (triangles in 2D) and by the parameters 
of the polynomial functions defined on the finite elements. The focus of the benchmarking is to compare the 
convergence speed rather than absolute runtime, i.e., the accuracy gain obtained by an increase in runtime. 
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Table 2. Computation times, real and imaginary parts and magnitudes of the 0*'' Fourier coefficients Ay{}tFC ~ 0) for 
polarizations TE aud TM and for the three bcriclimarkcd methods. Geometrical and material parameters arc denoted in 
Table 1 (data set 4). For each method, increased computation time corresponds to a higher spatial resolution. Please 
note that with Delight, both, TE and TM modes are computed simultaneously, JCMharmony was used with regular and 
with adaptive grid refinement. Units of the Fourier coefficients are [V/m] (TE), resp. [A/m] (TM). 



Therefore, other settings of the simulators were conservatively chosen in order to avoid any influence on the 
accuracy. This means, we expect that a more aggressive tuning of these parameters or implementations might 
slightly reduce the absolute runtime for all methods. 

As output we monitor the coefficients of the Fourier decomposition of the solution at the output boundary of 
the computational domain. The square of the Fourier coefficients is proportional to the power of light diflfracted 
into the corresponding diffraction order of the periodic mask. 



The Fourier coefficient of the y-component of the investigated field f = E, resp. H is defined as 



Ayikpc) = — fv{x,y,zo)exp{-ikFcx)dx, (5) 

Px J-p^/2 

where kpc is the projection of the wavevector of the investigated diffraction order onto the x—y— plane {kpc = 
for zero order and perpendicular incidence). Table 2 lists the Fourier coefficients for TE and TM polarization, 
obtained from the solutions of the the three bcnchmarkcd solvers with different resolutions. JCMharmony has 
been used in adaptive and in regular grid refinement mode (see Section 3.3). To ease the use of Table 2, we have 
marked the already converged digits in bold. 
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Figure 6. Relative error (AAq) of the 0*'' complex Fourier coefficient vs. normalized computation time for TM-polarization 

(a) and TE-polarization (b). Circles correspond to results obtained with JCMharmony (with adaptive (a), resp. regular 

(b) grid refinement mode), triangles to Delight and squares to SOLID E. The corresponding data is listed in Table 2. 

Figure 6 shows the convergence of the different methods. Plotted is the relative error of the 0*'' order complex 
Fourier coefficient, Ay{0), of the simulated field components, 



AAn = 



\Ay{0) - Ay{0)\ 

1^(0)1 



(G) 



where Ay{0) denotes the complex Fourier coefficient computed at highest spatial resolution with the correspond- 
ing method, vs. computation time. Since the computations were carried out on different platforms, computation 
time is given in units of a Matlab FFT run on the same platform. One unit of time corresponds to approx. 0.25 sec 
on a personal computer (Intel Pentium IV, 2.5 GHz). In spite of the respective internal convergence, we observe 
that the intensities of the 0*'' order diffraction computed with the three methods differ significantly (> 2%), 
especially for TM polarization. For that reason, we benchmarked the FEM solver using an example providing 
an analytical solution (see Section 5). The speed of convergence of the three methods differs also significantly. 

It can be seen from Fig. 6 (or, alternatively, Table 1) that the convergence behavior of the benchmarked 
solvers differs for TE and TM polarization. In the case of TM polarization, the FEM solver JCMharmony is 
the only one to reach target accuracies AAq < 10~^. For a low target accuracy of A^o > 10~^ computation 
times arc comparable for all three solvers, for a moderate target accurracy of AAq « 2 • lO^'' JCMharmony 
is about 100 X faster than Delight. SOLID E reaches an accuracy of about A^o ~ 10"'^. For low resolutions 
(computation time < 10 units) the loading of the code requires a significant amount of the total computation 
time (for SOLID E also the loading of the GUI which takes about 6 time units). In the case of TE polarization 
the convergence behavior of Delight is better than in the TM case, being slightly faster than JCMharmony at low 
target accuracies AAq > 1%, comparable at intermediate target accuracies AAq » 1%, and about 20 x slower 



than JCMharmony at high target accuracies of A^o ~ 10^^. Please; note that in the TM case adaptive; grid 
refinement yields a gain in convergence for the FEM solver JCMharmony, while in the TE case it is not superior 
to regular grid refinement (see also Section 3.3). 

5. BENCHMARKING WITH AN ANALYTICAL SOLUTION 

The attainable absolute accuracy of the finite element solver was assessed using an example which is close 
enough to lithography applications and allows for a fully vectorial 2D quasi-analytic solution by means of a 
series expansion with proven convergence. The geometry is depicted in Fig. 7. It consists of dense lines and 
spaces in an infinitely thin membrane mask of perfectly conducting metal embedded in free space. The metal 
layer covers the a:y-plane, the strips of the grating are oriented in y-direction. The grating is illuminated from 
the top by a TM polarized plane wave, i.e., Hy^ia = exp(— ifczZ — itot). Geometry plus incident wave imply 
the following properties: (A) Hy does not depend on the y-direction which allows us to write Hy = Hy{x,z), 
Maxwell's equations separate for the iJ^-component yielding a Helmholtz equation, (B) the field is periodic in 
X with period 2L, (C) the condition Hy[x,0) = 1 holds true inside the gap, [.t| < a, whereas the condition 
dnHy{x,0) = holds true outside this gap, a < |x| < L. Altogether, the time harmonic Maxwell's equations for 
the magnetic field component Hy reduce to 

dxxu + dzzu + k'^Hy = for (.x, z) E [-L, L] x 

Hy{—L,z) — Hy{L,z) (periodicity in x) 
dnHy{x,0) = fora<|a;|<L 
Hy{x,0) = 1 for \x\ < a 
+ radiation condition for the scattered field. 

Periodicity in x justifies the Fourier expansion Hy{x, z) = J2-n^ Cn{z) exp{imrx/L). Inserted into (7) we obtain 




Figure 7. Left: Geometry of the problem reduced to a 2D scheme in the x, z plane. Right: Intensity plot of the 
superposition of incident and scattered field (black: low intensity, white: high intensity). 

the ordinary differential equation c^(^) = ((titt/L)^ — A;^)c„(z) for the unknown coefficients c„(z). It has the 
general solution 

I \ _\ a„ exp(A/ {nir/LY — k'^z) + bn exp(--\/(n7r/L)^ — k'^z) for P < {nn/L)^ 
" \ anexp{i^/k'^ — {nn/LYz) + 6„ exp(-i-\/fc2 - {nn/LYz) for (n-KjVf < kP' 



10" 3 4 5 

10 10 10 

Number of Unknowns 



Figure 8. Convergence of the FEM approach towards the analytic solution. (See original publication for images with 
higher resolution.) 



The coefficients a„, 6„ have to be chosen such that a proper radiation condition holds true. This imphes 6„ = 
since we firstly require decaying solutions in case of (riTr/L)^ > fc^ and z — > — oo, and secondly, outgoing 
solutions in case of {m:/L)'^ < fc^ and z — > — oo. Based on the ansatz for Hy, the derived explicit solutions for 
the coefficients c„ in terms of coefficients a^, an uniform sampling, a;„ = —L + nL/N, n = —N, . . . , N — 1, we 
obtain 2N equations for 2N unknowns a„: 

Hy{xn, 0) = 1 = ^ ane*"""""/^ for n such that |a;„| < a 

n 

dnHy{xn, 0) = = ^ anVTW^P^^e*"^^""/-^ for n such that a < |a;„| < i . 

n 

The numerical solution of this system supplies the quasi-analytic reference solution Hy{x,0) with \x\ < L. The 

reference solution was computed with N = 2^"^ ensuring an error less than 10~^ with respect to the mean 

quadratic error, ||-ffy,Ar - i?y,oxact|| ^/N\/j2n {Hy,N{xn) - Hy^cKiictixn)f , on the interval -L < y < L aX 
z = 0. The reference solution was then compared to FEM solutions with increasing number of unknowns, where 
exactly the geometry depicted in Fig. 7 was used. The convergence results arc presented in Fig. 8. The first 
FEM solution was obtained for 587 unknowns yielding an error 1.506 ■ 10^^ , the last solution was computed 
with 82 303 miknowns yielding an error 3.158 • 10^*. This proves that the FEM solver JCMharmony produces 
simulation results which converge to the exact solution of the scattering problem. The current implementations 
of the other solvers did not allow to examine this example with these. 



6. CONCLUSIONS 

We have benchmarked a FEM solver for mask simulations against two other rigorous methods (FDTD and 
waveguide method). The FEM solver allowed to reach high accuracies which were not accessible with the other 
methods. The computation time for phase mask simulations at moderate target accuracies was up to two orders 
of magnitude lower when using the FEM solver compared to the other methods. Further we have performed a 
benchmark of the FEM solver against an analytically accessible problem, and we have shown the wide range of 
applications of the solver by examining several typical simulation problems. 
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